Draft version May 6, 2013 

Preprint typeset using I^TeX style emulateapj v. 5/2/11 



PROTOSTELLAR DISK EVOLUTION OVER MILLION- YEAR TIMESCALES WITH A PRESCRIPTION FOR 

MAGNETIZED TURBULENCE 

Russell Landry 

Physics Department, University of Texas at Dallas 



Sarah E. Dodson-Robinson 

Astronomy Department, University of Texas at Austin 



Neal J. Turner 

Jet Propulsion Laboratory/California Institute of Technology 
AND 

Greg Abram 

Texas Advanced Computing Center, University of Texas at Austin 
Draft version May 6, 2013 

ABSTRACT 

Magnetorotational instability (MRI) is the most promising mechanism behind accretion in low-mass 
protostellar disks. Here we present the first analysis of the global structure and evolution of non-ideal 
MRI-driven T-Tauri disks on million-year timescales. We accomplish this in a 1+lD simulation 
by calculating magnetic diffusivities and utilizing turbulence activity criteria to determine thermal 
structure and accretion rate without resorting to a 3-D magnetohydrodynamical (MHD) simulation. 
Our major findings are as follows. First, even for modest surface densities of just a few times the 
minimum-mass solar nebula, the dead zone encompasses the giant planet-forming region, preserving 
any compositional gradients. Second, the surface density of the active layer is nearly constant in time 
at roughly 10 g cm~^, which we use to derive a simple prescription for viscous heating in MRI-active 
disks for those who wish to avoid detailed MHD computations. Furthermore, unlike a standard disk 
with constant-a viscosity, the disk midplane does not cool off over time, though the surface cools as 
the star evolves along the Hayashi track. Instead, the MRI may pile material in the dead zone, causing 
it to heat up over time. The ice line is firmly in the terrestrial planet-forming region throughout disk 
evolution and can move either inward or outward with time, depending on whether pileups form near 
the star. Finally, steady-state mass transport is an extremely poor description of fiow through an 
MRI-active disk, as we see both the turnaround in the accretion fiow required by conservation of 
angular momentum and peaks in M{R) bracketing each side of the dead zone. We caution that MRI 
activity is sensitive to many parameters, including stellar X-ray fiux, grain size, gas/small grain mass 
ratio and magnetic field strength, and we have not performed an exhaustive parameter study here. 
Our 1+lD model also does not include azimuthal information, which prevents us from modeling the 
effects of Rossby waves. 



1. INTRODUCTION 

A typical, low-mass T-Tauri star accretes mass at a 
rate of 1Q ~^ M^ yr~^ — one Jupit e r mass every 100,000 
years ( e.g. [Hartmann et a l.| p!998 Sicilia-Aguilar et al. 
|2Q04l ). [Balbu s ^ Hawley (iMTjlput forth the magne- 
torotational instability (MRI) as the most likel y driver^ 
of accretion 




physical description of angular momentum transport 
was available , accretion was often mode led using the a- 
prescription ( Shakura fc Syunyaev||1973 ), 

(1) 



ac.H. 



The a-prescription relates turbulent viscosity to length 
and velocity scales in the disk based on dimensional anal- 
ysis. In Equation [l] v is the turbulent viscosity, Cg is the 
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sound speed and H is the pressure scale height, a is a 
dimensionless efficiency that is often assumed, without 
physical motivation, to be uniform throughout the disk. 

The first numerical investigations of MRI-dri ven tur- 
bulence wer e local shearing box simulations (Hawley 
et al. 1995), which treat a box of approximate size 
TttTTxTTx H centered at a given distance from the 
star. The shearing box is still a useful technique for in- 
vestigating detailed properties of turbulence. However, 
more recent investigations computed turbulent viscos- 
ity from first principles for global disk models, which 
were previously the domain of the a-prescription. Such 
global simulations confirmed that the MRI can lead to 
turbulence-driven accretion matching obs erved rates of 
10"^ Mc7) yr~^ in either unstratified disks (| Hawley || 200l " 



Steinacker fc Papaloizo u 2002 ; Lyra et al. 2 008p or thin 
stratitied disks (Sorathia e t al.||2010| ). MRI turbulence 
is self-sustaining for the simulation ti meframe of around 
1000 orbits near the inner boundary (Fromang fc Nelson 



2 



2006} [MockeTal]!^^ 

In parallel with the development of global accretion 
disk models came investigations of the behavi or of the 
MRI in non-ideal, partially ionized fluids. I Gammie] 
([1996^ was the first to point out that protostellar disks 
nave surface layers ionized well enough to couple to mag- 
netic fields, and an interior dead zone where extremely 
low ionization levels prevent magnetically driven turbu- 
lence. Subsequent investigations incorporating Ohmic 
resistivity into the MHD equations confirmed Gammie's 
prediction that a high enough^ resistivity would dampen 
the growth of the MRI ( |Jin[l996| |Sano & Miyama 1999). 
Numerical simulations showed a "dead zone" near the 
midplane and near the star, where the surface density 
is high en ough to shield the di s k interior from io nizing 
radiation ([Flemi ng et al.||2QQQ[ |Sano et al.||2QQQ|). Re- 
cent research has revealed that the dead zone is "not en- 
tirely dead: the shielded interior still experiences a shear 
stress only about an order of magnitude less than the 



activ e layer due to propagating acoustic waves (Fleming 
& Stone 2003| and smooth, large-sc ale magnetic fields 
fTOmer et ar[|2QQ7[ [Turner fc Sano||2QQ8D . Other non- 
ideal ettects that attect the MKi include the Hall current 

Sano & Stone 



Wardle 1999 



2002T 



Balbus < 



ar dittusion 



and ambi po 
Mac Low et al.|[T995l |Kunz 



Terque m 2001 

e^Blaes & Balbus.1994 



Balbus||2 004; Ba i & '^tonej 



201 ip (See Section 2.1 for a more complete description of 
each effect). Improved understanding of both non-id eal 
MRI and ionization in disks (Igea & Glassgold 1999 [iT 



gner fc Nelson^^2006, ) has provided the tools to descrioe 



the MKI in gas of any ionization fraction or density that 
can be found in a protostellar disk. 

Detailed short-timescale snapshots have now been con- 
structed of angular momentum transport in a protostel- 
lar disk , including the dead zone, turbulent layers and 
corona ( Dzyurkevich et al.|2010 Kretke & Lin|2010 Baij 



2011[ jFlaig et 
active I - Ihuri 



disks may 



Yet it is possible tnat IVlKI- 



may not be in steady state due to 
(a) the dead zone and its resulting, radially varying ac- 
cretion rate and (b) the lack of a protostellar envelope to 
provide material to maintain steady, inward mass trans- 
port. Our protostellar disk snapshot, then, must change 
over the million-year timescales on which the star/disk 
system evolves. The goal of this work is to illustrate how 
disk interior structure and angular momentum transport 
change over the entire, multi-million-ye ar lifetime of the 
T-Tau ri p hase. We extend the work of [Armitage et al^ 
( 2001 ) and Zhu et al. ( 2010 ) , who also performed million- 
year simulations, by calculating a radially and vertically 
varying a based on the ionization state rather t han as- 
suming a constant a- value in the active zone. 'Martin] 
|et al. (2012b) also modeled FU Orionis outbursts using 
time-dependant global simulations of MRI- active disks, 
including Ohmic resistivity, using 1-D layered models 
which used the a-prescription for the active layer, and 
an analyt ical approximation for the active layer surface 
density (Martin et al.|2012a[). Our models build on their 
work by adding a vertical dimension to the disk struc- 
ture. 

Our evolving model of a magnetically turbulent T- 
Tauri disk answers the following questions: 

1. How do the relative sizes of the dead zone and ac- 
tive layers change over time? 



2. How does M vary with radius and time? 

3. How can disk modelers parameterize heating in the 
active layers and dead zone without resorting to a 
3-D MHD simulation? 

4. Does the disk midplane heat up or cool off over 
time? 

5. Where is the ice line in an MRI- active disk and how 
does its location change over time? 

Questions 1, 2 and 3 elucidate basic properties of an 
MRI-turbulent accretion disk. Questions 4 and 5 high- 
light fundamental ways in which disk models based on 
the standard viscosity prescription with constant a lead 
us astray. Note that our disk model does not include 
photoevaporation, which is another process that oper- 
ates on million- year timescales that can produce radially 
varying accretion accretion rates. 

Our ability to simulate an entire T-Tauri disk lifetime 
is due to a new MRI activity prescription that allows us 
to compute the thermal and viscous effects of MRI tur- 
bulence without resorting to 3-D magnet ohydro dynamic 
simulations of the turbulence itself. We can thus reduce 
our computational domain from three spatial dimensions 
to 1+1 spatial dimensions — 1-D vertical structures rep- 
resenting axisymmetric disk annuli that are connected 
only by a 1-D radial m ass transport equation (Dodson- 
Robinson et al. 2009). Sacrificing information about 
small-scale turbulent liuctuations, we retain our ability 
to accurately describe large-scale structures such as the 
dead zone and active layers while dramatically improving 
our ability to simulate long timescales. 

Our paper is organized as follows. In Section [2] we dis- 
cuss the basic equations governing the MRI under non- 
ideal MHD conditions and give our prescription for de- 
termining turbulent viscosity. We present our method 
of computing vertical strtucture and mass transport in 
Section |3] We outline a basic picture of MRI-turbulent 
disk evolution and answer Questions 1, 2 and 3 in Sec- 
tion |4j In Section [5) we discuss the differences in thermal 
structure between our model and constant-a disk models, 
which leads us to answer Questions 4 and 5. We discuss 
the limitations of our model in Section |6] and present 
our conclusions in Section [7j Readers who wish to skip 
over the details of the computations may wish to proceed 
directly to Section [Ij 

2. SIMULATING MAGNETOROTATIONAL 
INSTABILITY-DRIVEN TURBULENCE IN PARTIALLY 
IONIZED GASES 

When an accretion disk is fully ionized and the mag- 
netic field is weak, the entire disk is MRI turbulent. Yet 
when a disk is only partially ionized, as is the case for a 
protoplanetary disk, there is an incomplete coupling be- 
tween the disk gas and the magnetic field, and non-ideal 
effects become important to the growth of MRI-driven 
turbulence. In this section, we describe how we treat 
angular momentum transport from MRI turbulence in 
non-ideal, partially ionized gases. We begin by describ- 
ing our turbulence criterion in Section [2?H then list our 
method for determining the di ffusi on regime and result- 
ing turbulent stress in Section [2^ 



2.1. MRI activity criteria in the three non-ideal regimes 

The non-ideal magnetic induction equation has three 
extra terms corresponding to the three non-ideal effects 
( [War die 2007^ , Ohmic resistivity, the Hall effect and am- 
bipolar dittusion: 



dB 

~dt 



V X X B) 



V X 



T^oV X B + r]H{^ X B) X B^ r]A{V x B)_ 



(2) 



In Equation l2j v is the gas velocity, B and B are the 
magnetic field and magnetic field unit vector, and ± 
refers to the component perpendicular to B. r]o^ Vh^ 
and rjA are the Ohmic, Hall, and ambipolar diffusivi- 
ties respectively. Ohmic resistivity dominates other non- 
ideal effects when collisions with neutrals cause both elec- 
trons and ions to decouple from field lines. When colli- 
sional drag is sufficient to decouple ions and grains from 
the magnetic field, but not electrons, the relative veloc- 
ity between the ions and electrons is non-negligible and 
the Hall effect is dominant. In the ambipolar diffusion 
regime, electrons and ions decouple from the neutral gas 
and the magnetic field lines are frozen to the charged 
species and drift through the neutral gas. 

When Ohmic resistivity is the largest non-ideal effect, 
the MRI will only occur if the Elsasser number. 



A 



(3) 



Stone 2002 



Turner 



is at least of order unity (iSano 

et al. 2007). In Equation ([3|, vaz is the Altveii speed 
in the vertical direction anot] is the Keplerian angular 
velocity. Physically, A is the ratio of the wavelength of 
maximum growth to the diffusive scale length. The tan- 
gled magnetic fields in MRI turbulence usually have a 
toroidal component with pressure 10 to 30 ti mes greater 
than the p ressure in the vertical component ([Turner fc 
Sano||2008[) , so the Alfven speed in the vertical direction 
used to calculate A is 



1 

10 



(4) 



where va = B / ^/Aifp is the total Alfven speed. The A > 
1 criterion ensures that the most unstable mode can grow 
more quickly than the charged parti cles can diffuse across 
magnetic field lines. Simulations by Sano & Stone (|2002 ) 
suggest the Hall effect does not change the conditions 
for turbulence if Ohmic diffusion is also present. Further 
work is required to determine the growth of the MRI in 
regime s where the H all term is much stronger than other 
terms ( I Wardle|[2QT2 ) . 

Ambipolar dittusion arises from the relative motion of 
ions and neutral particles in the disk gas. In the "strong 
coupling" limit, in which ion density is negligible and 
electron recombination time is much smaller than the 
orbital tim e as is generally the case for protoplane- 
tary disks ( Bai|2011 ), ion density cannot be assumed to 
follow the continuity equation. Instead, the ion density 
is determined by the ionization-recombination equilib- 
rium, and characterized by the parameter Am (Chiang 



Murray-Clay|[2QQ7| ): 



Am = 



(5) 



where pi is the ion density and 7 is the neutral-ion drag 
coefficient, 

1=- r-I-' (6) 



rur, 



In Equation [6) (Jni is the effective cross section for 
neutral-ion collision and Wni is the relative velocity be- 
tween neturals and ions. Physically, Am is the ratio of 
the orbital period to the collisional times cale between 
ions and neutrals. Since rjA = v\/^pi ( |Bai fc Stone 
2011), one can rewrite Equation [s] as 



Am 



(7) 



which is equivalent to the Elsasser number A in the 
Ohmic regime. Similarly then. Am is the ratio of the 
wavelength of maximum growth to the ambipolar diffu- 
sive scale length. 

In their three-dimensional shearing-box simulations ex- 
plorin g the effect of ambip olar diffusion on MRI turbu- 
lence, Bai & Stone (2011) determine that heavily ion- 
ized yet tenuous disks can only sustain turbulence when 
threaded by weak magnetic fields. The magnetic field 
strength is characterized by the plasma /3, the ratio of 
the gas pressure to the magnetic pressure: 



/3- 



SirP 



(8) 



The requirement that the magnetic field energy be small 
in comparison to the gas thermal energy (a "weak" field) 
restricts MRI turbulence to v alues o f the plasma /3 that 
are greater than a minimum (Bai & Stone 2011): 



/3min{Arn) 



50 



Am^ 



Am^-^ 



1/2 



(9) 



The maximum field strength, beyond which the field is 
too strong to be destabilized for any given field geom- 
etry, decreases the more important ambipolar diffusion 
becomes (smaller Am). However, for a sufficiently weak 
field, MRI can be sustained even for Am <C 1. 

In a protoplanetary disk, ambipolar diffusion domi- 
nates in the atmosphere, which is diffuse and highly ion- 
ized by stellar X-rays and cosmic rays. Ambipolar dif- 
fusion is also important in the outer disk where the sur- 
face density is very low. Ohmic resistivity dominates in 
the dense inner region shielded from ionizing radiation. 
Ignoring the effects of the Hall diffusivity, which are un- 
likely to alter either the conditions required for MRI or 
the strength of the turbulence wh ere Ohmic dissipation is 
also present (IS ano fc Stone|20'Q2 ), a non-ideal protoplan- 
etary disk differs from an ideal accretion disk through 
the possibility of a dead zone in the inner disk. The 
dead zone would remain cold and would not efficiently 
transport angular momentum. The accretion efficiency 
in the upper, ionized layers of non-ideal disks also lags 
behind ideal disks due to ambipolar diffusion effects. 

In order to compute the disk viscosity and angular 
momentum transport properties, we need to know (a) 



4 



whether the MRI is operating, and (b), if so, how strong 
the turbulence is. Since the MRI growth timescale is 
roughly 1/1], hundreds of thousands of times shorter than 
our ~ 1 Myr simulation timescale, we assume the MRI is 
either fully saturated or completely damped. Our MRI 
turbulence criterion is therefore equivalent to that of|Bai] 
(I20TTI): 



actions: 



• If A > 1 and /3 > Prnim neutral gas couples to the 
magnetic field and MRI is saturated; 

• If A < 1 or /3 < Prnin^ ucutral gas decouples from 
the magnetic field and MRI is damped. 

In the next section, we discuss the computation of the 
magnetic diffusivities that determine A and Am. 

2.2. The diffusion regime and turbulent stress 

We cannot apply our turbulence criterion without 
knowing the values of 7^0 , Vh and ija- For a given charged 
species j, the ratio of Lorentz force to the neutral drag 
force is 



(10) 



where Zje is the charge of j (negative or positive), B 
is the magnitude of the magnetic field, mj is the mass 
of j, c is the speed of light, and 7^ is defined according 
to Equation Im (Note that Pj is not the same as the 
plasma (3 of Equation [s]) For each diffusion regime, one 
can define a conductivity by summing over all charged 
species: 



CFp 



ec s-^ ^3^3 



B 



^ 1 



(11) 

(12) 
(13) 

(Wardle 2007). In Equations 11-13, Uj is the number 
density ot species j. Finally, one can write the diffusivi- 
ties according to 



B 



^3 ^ 3 ^3 

PI 



Vo 



VA 



4:7T(Jo ' 
^2 



^2 



47r<Jx cr± 



(14) 
(15) 
(16) 



where aj_ 



/p. 



To determine the equilibrium abundances of charged 
species nj , we s olve a simplified set of c hemical reactions 
from Model 4 of II gner & Nelson rt2006fc , which we briefly 
motivate here. The set is derived trorn the following re- 



H2 + X ^ H+ + e- 



H 



H2 



H+ + Ha 

H+ + CO ^ HCO+ 
2H + g ^ Ha + g 
HCO+ + e- ^ CO + H 
HCO+ + Mg ^ Mg+ + CO 

Mg+ + e" -> Mg, 



H 



(17) 
(18) 
(19) 
(20) 
(21) 
(22) 
(23) 



where HCO+ is a representative molecular ion, Mg+ is 
a representative metal ion and g is a grain. Here every 
species (except the energetic particle X — a cosmic ray, X- 
ray or radionuclide decay product) is created in at least 
one reaction, and destroyed in at least one other. Over 
the whole set, no species is produced or consumed on 
balance. The subset producing the ions and electrons 
reduces to 



2H2 + 2X + 2C0 ^ H2 + 2HC0^ 



2e- 



(24) 



That is, each energetic particle striking a hydrogen 
molecule yields one ion and one electron. In constructing 
the conductivity lookup tables we therefore approximate 
Equations 17-21 by 



H2 + X ^ HCO+ 



HCO^ 



(25) 
(26) 



neglecting the fact that the molecular ion contains just 
one hydrogen atom. Since HCO+ is orders of magnitude 
less abundant than H2 , forming ions leaves the H2 density 
unchanged. Similarly, we don't model CO destruction 
and reformation because the ion is so much less abundant 
than the molecule. Equation 22 then becomes 



HCO+ + Mg ^ Mg+ + H2 



(27) 



The simplified network consists of Eqs. 25, 26, 27 and 
23, together with the grain surface reactions described by 
Ilgner & Nelson (2006). The metal atoms' ther mal ad- 
sorption and des orption on the grains is included. |Ilgner| 
Nelson (2006) found that this reduced network yields 
similar results to a detailed version including hundreds of 
species and thousands of reactions, in the most common 
situation where the recombination occurs mostly on the 
grains. 

The internal grain density, gas/small grain ratio, and 
grain size used in the chemical reaction network are listed 
m TablelH Here we deviate from the standard interstellar 
gas/small grain mass ratio of 100 and assume some grain 
growth has occurred, so that 90% of the grain mass is 
in grains larger than 1/im. Using the standard gas/dust 
ratio of 100 resulted in no MRI turbulence (see section 
4.1 ). To avoid having to run the chemical reaction net- 
work at every timestep of o ur million- year sim ulations, 
we followed the approach of |Flaig et al. (2012) and cre- 
ated a look-up table of magnetic ditfusivities as a func- 
tion of temperature T, gas density p, ionization rate ( 
and plasma p. To generate the look-up table, we ran 
the chemical reaction network until it reached equilib- 
rium abundances of all species for each combination of 
T, p, C and /3. We then computed the conductivities and 
tabulated diffusivities according to Equations 10-16. 
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The rate coefficient for tiie reaction H2 + X — ^ H2" + e~ 
is, of course, the ionization rate (. To determine we 
consider cosmic rays, stellar X-rays, and short-lived ra- 
dionuclides (ri/2 <C 10^ yr), including ^^Al. Follow- 
ing Ume bayashi fc Nakano] (|2009), we take the ioniza- 
tion rate from short-lived radionuclides as (r = 7.6 x 
10~^^ s~^ and calculate the attentuated cosmic ray ion- 
ization as: 



CcRiz) 



^surf 
^CR 



exp 



-\-exp 



^CR 



-3/4" 



-3/4 



4/3 ^ 



(28) 



where Cci?*^ unattenuated cosmic ray ioni zation 



rate, Xqr is the co smic ray penetration depth (Ume- 



bayashi fc Nakano [l981), and 111^2 (^) are the mass 
columns above and below the vertical height z. Values of 
C(jRi ^CR and all other numerical inp uts to our model 
are listed in Table [l] Finally, following "Bai & Goodman 
( [2009 J , we calculate the stellar X-ray ionization rate Cx 
as: 



Cxiz) 



(Ml 



(29) ■ 



where Lx,29 = ^x/(10 erg s~ ), and Lx is the the 
stellar X-ray luminosity. We take Lx,29 = 20 to match 
the young solar-mass stars observed in the Orion Nebula 
( Garmire et al.|200 0V Here we keep the stellar X-ray flux 
constant m time, though it could certainly vary in either 
a smooth, systematic way with age or stochastically with 
accretion bursts. All parameters in Equation|29]are listed 
in Table [1] The first exponential represents attenuation 
of X-rays by absorption, and the second represents the 
contribution from scattered X-rays. We show vertical 
profiles of the ionization rate at two different disk radii 
in Figure [1] 

After the ionization rate is determined at each (i?, z) 
zone in our disk model, an interpolation through the 
look-up table can return the magnetic diffusivities. We 
can then compute the Elsasser number A and Am and de- 
termine whether MRI turbulence in the zone is active or 



not, according to our turbulence criterion (Section 2.1). 
The last remaining ingredient in our MRI prescription 
is a rule for determining the strength of the turbulence, 
where it is present. Where MRI is saturated, we use the 
scaling relations 



1 

ac.H ■■ 



(30) 
(31) 



found between turbulent stress and magnet ic field 
strength in a variety of s hearing-box simulations (Hawley" 
|et al.„199H |Sano et al.|[2Q04, Bai fc Stone, 2011} . Note 



-4/3 



that Equation 31 applies no matter the field geometry or 
value of A or Am. 

Since there is an accretion flow cause d by large-scale 
magnetic fi elds even in the d ead zone (iTurner fc Sano 
|2008[ jNelson fc Gressel 2010), we set a minimum value 
of a where the MRI is damped. In the active layer, a is 
close to its maximum value of ~ 0.5, set by the cessa tion 
of MRI in the strong- field limit o f /3 > 1 dHawley et 
lL| [l995l | B5bus fc Hawl^[T998 |. |Bai fc StoSil ^m\\ 



found a similar result in the ambipolar regime, a ^ 0.4 
for Am 00. The shear stress in the dead zone is an 



order of magn i tude less than in the active layer ([ F leming 
Stone|[2QQ3l [Turner et al^[2QQ7| [Turner fc Sano||20n5l r 
d the plasma b at the midplane is typically two t( 



and the plasma (J at the midplane is typically two to 
three orders of nia^nitude higher than at the top of our 
grid (see Figure ^ which shows profiles of the vertical 
component of /3, /^^ = SttP/^^, for two different disk 
radii), amin must be therefore be roughly four orders 
of magnitude less than a^ax to produce an appropriate 
level of dead zone shear stress. We take amin = 10~^ 
(Table [1]). 

3. DISK STRUCTURE AND MASS TRANSPORT 

With our method of computing turbulent viscosity at 
any location in a protostellar disk, we may now simu- 
late how the disk re-distributes its mass throughout its 
multi-million- year existence. Since our viscosity prescrip- 
+tion depends on the four inputs T, p, (3 and which are 
functions of radius and vertical height {R^z)^ we must 

} compute the detailed vertical and radial structure of the 
disk. Our computational setup, simila r to the disk mod- 



els of Dodson- Robinson et al. (2009), is based on the 



following simplifying assumptions: 

1. The disk is axisymmetric and symmetric about the 
midplane; 

2. The disk is geometrically thin, H/R <C 1; 

3. Heat escapes in the vertical direction much faster 
than it is carried with the gas flow in the radial 
direction. 

Assumption 1 reduces the physical three-dimensional 
disk to a two-dimensional quadrant with zero flux at the 
midplane (required for symmetry). Neglecting the az- 
imuthal dimension is a critical step in speeding up the 
code to allow long-timescale simulations. Assumption 2 
allows the vertical and radial dimensions to be decoupled 
in a 1+1-D framework, so that energy transport proceeds 
only in the vertical direction. Each radial gridpoint con- 
tains an independent vertical structure model. Assump- 
tion 2 is valid as long as the vertical sound-crossing time 
(of order the orbital timescale) is much less than the 
accretion timescale, which is generally true in T-Tauri 
disks. 



In section [3.1 [we des cribe our vertical structure model, 

dis cuss our mass transport 
3.3| contains a description of 



while in Section 
parametrization 



|3.2[ we 

Section 



our computational methods 

3.1. Vertical structure 

Hydrostatic equilibrium and the thermal balance be- 
tween stellar irradiation, radiative cooling and viscous 
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Fig. 1. — Plots showing the multi-layered structure of evolved, magnetically active disks. Besides the dead zone and the active layer, 
both disks have a corona near the star and at the surface where the magnetic field is too strong for the tenuous gas to create turbulence. 
Disks may also have a double-layered active zone with a dead slice in some locations. Top left: a(z), pz(z) and ((z) of Model 1 at 1.5 AU 
after 1 Myr. Top right: Model 1 at 10 AU. Bottom left: Model 2 at 1.5 AU, 1 Myr. Bottom right: Model 2 at 10 AU. In each panel, 
the active region with A > 1 and ^ > ^min is shaded with red dots. 

Vad = 2/7 is the adiabatic thermodynamic gradient for 
diatomic gas and 



heating govern the vertical structure of our disk. We use 
the flux-hmited diffusion approximation for the transport 
of viscously generated energy. The accretion flux gradi- 
ent is determine d by the visco us energy generation rate 
per unit volume ( Pringle||1981 ): 

dFacc 9 ^^^^ 



3 tzPF 
~iacQ?zT^' 



(38) 



dz 
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In Equation [32j Vt is the height-dependent Keplerian fre- 
quency 

1/2 

GM* 



is the radiative thermodynamic gradient, where a is the 
radiation density constant and n is the local Rosseland 
mean opacity. For full details on how to compute V conv 
in the disk's convective zone, see|Kippenhahn fc Weigert 



+ ^2 



,3/2 



(33) 



Temperature and pressure are related by the ideal gas 
equation of state, 



(34) 



where ji = 2.33 g mol~^ is the mean molar weight of 
the disk gas (Table [T]). The temperature and pressure 
gradients are 

dTacc dP T 



dP 

dz 



(36) 



In Equation [35j Tacc is the temperature contribution 
from viscous heating only — stellar irradiation and the 
ambient molecular cloud also contribute some of the ther- 
mal energy. T is the true temperature and includes all 
heat sources. 

To calculate V = dlnT/dlnP, the thermodynamic 
gradient, we use the Schwarzschild criterion for stabil- 
ity against convection: 



(1994). Computing V requires the local Rosseland mean 
opacity, tv. At low temperatures, T < 700 K, we use 
the opacities of |Semenov et al. (2003) calculated for a 5- 
layered sphere topology. At higher temperatures where 
molecular ^a s dominates opacity (T > 1000 K), we use 
the tables of [Ferguson et al.| ( |2005D . For 700 K < T < 
1000 K, we interpolate between the two tables using a 
weighted average in log(T) space. For more details, and 
plots of the resulting opacities, see |Dodson- Robinson et 
al. (2009). 

Integrating the coupled ODEs in Equations 32, 35 and 
36 requires two different temperatures: T, the true tem- 
perature resulting from all sources of thermal energy, and 
Tacc^ the component from viscous heating only. The 
other heat sources are the central star, which sets an 
equilibrium temperature component Tgg, and the ambi- 
ent star-forming region, from which long- wavelength ra- 
diation that penetrates the disk sets a minimum temper- 
ature Tamb' To compute Teg, we begin by assuming the 
disk surface is flared as a result of hydrostatic equilib- 
rium and radiative equil ibrium with the star. Following 
the models developed by Chiang & Goldreich ( 1997J , we 
calculate the grazing angle at which stellar energy en- 
ters the disk: 



8 /t; 
7 [n 



4/7 



2/7 



(39) 



^rad^ ^ rad ^ ^ ad 
^ conv •) ^rad 



(37) 



where i?* and T* are the star's radius and effective tem- 
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perature and Tc is a measure of the gravitational poten- 
tial at the surface of the star: 



crR^ 



(40) 



In Equation [40j is the star's mass and a is the 
Stefan-Boltzmann constant. As the disk evolves, we de- 
termine the star's temperature and radius as a function 
of age from the pre main - seque nce evolutionary tracks 
of D'Antona & Mazzitelli (1994). Since the star's lumi- 
nosity decreases as it moves down the Hayashi track, the 
disk flaring becomes less pronounced over time and the 
disk surface, where T ^ Tg^, cools. 

In the direction parallel to the disk midplane, the stel- 
lar radiation penetrates to an optical depth r| ~ 1. 
The asterisk subscript denotes that this optical depth 
is measured at the peak wavelength of the starlight, 
near 1 /im. Measured perpendicular to the disk mid- 
plane, stellar radiation is mostly attenuated by optical 

depth ^ tIo ^ 0. The equilibrium temperature be- 
tween the stellar heating and radiative cooling at the disk 
surfa ce — neglecting visc ous heating — is ( [P'Alessio et al.| 
'20061 iNatta et all 20001): 



^ eq,surf 



: 0.8 



^surf 



l/^ /-n \ 1/2 



(41) 



where Tsurf is the Rosseland mean optical depth at the 
disk surface for blackbody radiation at Teq^surf- The top 
surface of our vertical grid is defined by Tsurf = 0-2. 
Since half of the stellar radiation absorbed by grains at 
the disk surface is re-radiated into space, the equilibrium 
temperature with the star as a function of height in the 
disk — again, neglecting viscous heating — is 



Te,{zf = \Tt,^,^,fe-^% (42) 
where is the Rosseland mean optical depth to height 



-surf 



np dz. 



(43) 



Finally, the true temperature T simply the flux sum of 
the individual temperature components: 



-'- acc ~^ eq ~^ amh' 



(44) 



3.2. Radial 



usion 



Since the radial and vertical dimensions of our disk 
model are decoupled, we must treat mass transport as a 
one-dimensional problem. Yet the MRI-active, partially 
ionized disk is vertically layered, with the most active 
accretion occurring at the surface. The key to a suc- 
cessful 1-D description of layered accretion is the fact 
that vertical re-distribution of mass within an annulus 
occurs more quickly than accretion in the active layers: 
l/r^ <C /inactive- By Computing a mass- weighted, ver- 
tically averaged value of turbulent viscosity in each an- 
nulus, 



2 r-^f 



vpdz^ 



(45) 



where S is the surface density in the annulus, we can de- 
scribe mass transport using the radial diffusion equation: 



dt 



RdR 



(46) 



In each (i?, z) zone, we compute viscosity according to 
Equation [l] We use a height-dependent modified scale 
height 

(47) 



H = 



Vl + (2^2^7c2)' 



softened into a non-singular form ( Milsom et al.l |1994 
In MRI-active zones, a is given by Equation |3l| whi 
in inactive zones we set a to our chosen value of amin 
(Table [T]). The sound speed used to compute u{R,z) is 

o R 



'-T. 



(48) 



3.3. Computational methods and initial conditions 

To initialize the disk evolution model, we compute the 
vertical structure of a disk with the following features: 

1. A surface density profile S ex predicted by 



Zhu et al. (2010) for layered accretion disks in the 
l-l'auri phase. 

2. A pre-main-sequence star with a mass of 0.95 Mq 
and an initial age of 0.1 Myr, which roughly co- 
incides with the beginning o f the T-Tauri phase 



(Dunham & Vorobyov 



tmue to accrete a sma. 



2012|. The star wih con- 
amount of mass from the 
disk during the ~ 3 Myr T-Tauri phase. 

3. An ambient temperature of 20 K (Table [Tl to 
match the typical background tempera tures 
frared dark clouds ( Peretto et al.|[2QTQ| ). 



4. An outer radius of 70 A U, set by the 80 AU 
solar nebula size limit of Kretke et al. (2012). 
(Note that the disk expands from its initial ra- 
dius.) Kretke et al. show that a solar nebula with 
Rout > 80 AU would excite Kozai oscillations in 
some of the planetesimals scattered by Jupiter and 
Saturn, stranding them in stable, high-inclination, 
low-eccentricity orbits that surveys have not de- 
tected. 



5. An inner radius of Rin = 0.5 AU. The requirement 
that the disk stay below the dissociation temper- 
ature of H2, so that the ideal gas equation holds, 
dictates our choice of Ri^. |Ruden fc Lin| ( |1986[ ) 

' Kin does not anect the 



6. 



find that the exact value of Rin does not 
overall disk structure as long as Rout ^ Rin- 

A disk mass of either O.O15M0 (Model 1) or 
O.O3M0 (Model 2). Our disk masses are designed 



ulations (e.g. 


Lyra et al. 


2008 


2010 


Bai|201 


l|), almost a 


il of w 



[Dzyurkevich et al. 
nich use minimum- 



are probably n ot viable giant planet- forming envi- 
ronments (e.g. Thommes et al. 2008). In a forth- 
coming study, we will examine the evolution of 
MRI-active, high- mass, planet-forming disks. 
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We then evolve the disk forward in time using Equation 
46} We use fully implicit finite differencing, adjusting the 
timestep At so that surface density varies by a maximum 
of 1.0% during a single timstep. The inner boundary Rin 
experiences zero stress, such that matter falls directly 
from Rin onto the star. The disk is allowed to expand 
freely from the outer boundary Rout , with four new zones 
added to the disk each diffusion time of Rl^^/u{Rout)' 

At each timestep, we independently calculate the ver- 
tical structure for each zone in the radial grid. We be- 
gin the vertical structure solver with initial guesses of 
Tacc^ P and vertical component of the magnetic pres- 
sure Pb,z = ^z/^^ a^ of our grid, defined by 
'^surf = 0.2. We find the height of the grid surface by 



^surf 



K{p,T)P 



surf 



(49) 



The accretion flux at the grid surface is Facd^surf) = 
(jT^cc- We use a fourth-order Runge-Kutta integrator 
with adaptive stepsize control to integrate the coupled 
ODEs in Equations 32j 35 and 36 from the surface to the 
midplane. The vertical magnetic pressure stays constant 
in height, though it varies with radius. A solved vertical 
structure model has the properties 



Facc{z = 0) 



0, 



(50) 



required to keep the disk symmetric about the midplane, 

^surf 



2 1 pdz = S, 

lz=0 



(51) 



SO that the volume densities add up to the surface density 
in the annulus, and 



P{z = 0) = 1000. 



(52) 



Equation 52 requires that the plasma (3 be constant 
throughout the disk midplane. After turbulence is es- 
tablished, Fromang fc Nelson^ (2006) find midplane val- 
ues of 100 ^ P ^ 1000 lor a wme range of vertical box 
sizes, resolutions and boundary conditions in global ideal 
MHD calculations. We use the Newton- Raphson algo- 
rithm ( Press et al.||1992 ) to adjust the initial guesses of 
Tacci P and Fb z until a solution is found that satisfies 
Equations p)|5ll and|52l 

The on-off nature ofthe MRI creates discontinuities 
in and u{z) and Facc(^) that can cause the Newton- 
Raphson algorithm to oscillate between two sets of input 
parameters that bracket the correct solution of Equations 



50, 51 and 52 To avoid such oscillations in the Ohmic 



regime, we decrease the value of a gradually in the range 
1.6 > A > 0.4, using a sigmoid function: 



«/(i + e-'=(A-i)), 



(53) 



where k = 1.5\n{l/ {a minP) — !)• We use a similar sig- 
moid function to smooth a in the ambipolar regime for 

<I3< l.Qprmn- 
4. MASS TRANSPORT IN MRI-ACTIVE DISKS 

Here we present our simulations of the evolu- 
tion of magnet icall y turbulent disks over million-year 
timescales. In §4.1[ we discuss the relative sizes of the 
dea d zo ne and active layers (Question 1 in Introduction). 
In §4.2[ we demonstrate how the shrinking of the active 



layer over regions of high density enhances mass pileup 
in the dead zone. We also analyze the radial mass flow 
M{R) through the disk and show that the disk never 
reaches a steady state, even on millio n-year timescales 
(Question 2 in Introduction). In ^4.3, we give a simple 
prescription for accretional heating due to MRI for use 
with semi- analytical, non-evolving disk models to predict 
observables (Question 3 in Introduction). 

4.1. Turbulent Structure 

Figure [2| shows viscosity as a function of (i?, z) for the 
inner 30 of Model 1 (O.O15M0, top) and Model 2 
(O.O3M0, bottom). Outside of 30 AU, viscosity is almost 
independent of height z. The left-hand panels of Figure 
[2] show the disks after 10^ years of simulation time — 
a star age of 0.11 Myr, since we began the simulations 
with a st ar at the beginning of the T -Tauri phase, age 
0.1 Myr ( [Dunham fc Vorobyov||2012D . The right-hand 
panels show both model disks alter 1 Myr of simulation 
time. The plots reveal two important features of the 
evolution of low-mass, MRI- active disks: 

1. A midplane dead zone, where Ohmic diffusion 
quenches the MRI and restricts viscosity, extends 
to 16 AU in Model 1 and 21 AU in Model 2. 

2. Although the vertical heights of both the dead zone 
and the overall disk shrink with time as the disk 
loses mass and cools (see Section [5|, the radial ex- 
tent of the dead zone stays approximately constant 
in time. 

The radial size of the dead zone is larger than the ~ 
5 AU typically quoted for disks simi lar to the minimum- 
mass solar nebula (MMSN) (Matsumur a fc Pu dritz|2003 

1 n -WTT n I I ^ ^ ^ ^ I < I . , I I -J / " Ml I I ^ 



Salmeron fc Wardle 2008 I^ laig et al.'2{M). Here we 



detine the dead zone as the region ol the disk where A < 
1 . The primary reason our dead zone is so extensive is the 
surface density of our disks: the MMSN contains O.OIMq 
within 100 AU of the sun, while our least massive disk 
contains O.O15M0 within 70 AU. Our dead zone is also 
extensive becaus e we treat ambipol ar diffusion, as well as 



Ohmic diffusion. Sano et al. (2000) use the same tenfold 



dust depletion we do, but treat only Ohmic diffusion and 
find a dead zone that extends 5-10 AU. 

Dead zones with a 10~^ are required for maintaining 
compositional gradients such as the ice gradient in the 
asteroid belt, which would radially diffus e on million- year 
timescales for fully turbulent disks (Nelson fc Gressel 



20101). The fact that Uranus and Neptune have atmo- 
spheres with higher carbon enrichment than Jupiter and 
Saturn provides some evidence that the solar nebula had 
compositional gradients covering the entire giant planet- 
formation region (Encrenaz 2005). Given that the dead 
zone in a disk of just O.Olsli© can reach ~ 15 AU, the 
outer boundary of the giant planet-fo rmation region in 
the Nice model (Tsiganis et al. 2005), the solar nebula 
could have supported such a comipositional gradient. The 
lack of change in the dead zone radius over time suggests 
that compositional gradients are stable over million-year 
timescales. 

Figure [l] gives further insight into the vertical struc- 
tures of Model 1 (top) and Model 2 (bottom). In each 
panel, the active layer is shaded with red dots. At 10 AU, 
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Fig. 2. — Viscosity as a function of (R, z) at two timepoints during disk evolution. The more massive Model 2 disk has the larger dead 
zone, which extends past 20 AU. In both models, the radial extent of the dead zone stays approximately constant with time. Top left: 
Model 1 (disk mass O.OISM©), 10^ years of simulation time (star age 0.11 Myr). Top right: Model 1 after 1 Myr of simulation time. 
Bottom left: Model 2 after 10^ years of simulation time. Bottom right: Model 2 after 1 Myr of simulation time. 



both disks show the classical layered structure of an MRI- 
active zone on tope of a dead zone. At 1.5 AU, however, 
we see additional complexity in the disk's vertical struc- 
ture. At 45 mG (Model 1) and 70 mG (Model 2), the 
magnetic tension is strong enough to prevent MRI from 
bending the field lines in the low-density surface gas de- 
spite a high ioni zation rate of 1Q~^^ s ~^, forming a sta- 
ble corona (e.g. jMiller & Stone 2000|. Our stabilizing 
magnetic field values are roughly consistent with those 
of Salmer on fc War die (^008), who calculated that un- 
stable MRl modes can only grow for B < SO mG in the 
presence of Ijam grains. (See ^2.2| and Table IT] for more 
on the grain properties used in our models.) Figure [s] 
is a zoom-in on the stress coefficient a{R^z)^ viscosity 
iy{R,z) and ionization rate ({R^z) in the inner 4 AU of 
Model 2 after 1 Myr. The inactive corona is apparent 
for radii R ^ 3 AU, but moves above our computed disk 

surface (the location where 1) as the magnetic field 
stength decreases with radius. 

Also noticeable in Figures [l] and [3] is a split active layer 
in the inner part of Model 2. MRI- active regions of high 
a and u sandwich a "dead slice" that has reduced stress 
and viscosity by an order of magnitude. (The fact that 
a does not immediately plunge to its minimum value 
in the dead slice is a result of the sigmoid smoothing 
described in ^3.3[ ) While not present at the start of our 
simulations, the split active layer appears after only 5000 
years of disk evolution and persists until the end of the 
simulation at 3 Myr. The split active layer extends from 
roughly 1.2 A\J < R < 1.7 AU (Figure [3|, though its 



radial extent shrinks slightly as the disk evolves. 

The dead slice, in the part of the disk where ambipolar 
diffusion is the strongest non-ideal MHD term, is the re- 
sult of two competing effects. First, MRI in the ambipo- 
lar regime requires high ionization: decreasing ( toward 
the disk midplane reduces Am and shuts down turbu- 
lence. Yet ambipolar diffusion is quenched at high den- 
sities: increasing p toward the midplane increases Am, 
favoring turbulence. In the dead slice, the dropping ( 
temporarily dominates over the increasing p and shuts 
down the MRI. Model 1 does not have a dead slice be- 
cause its surface density is about 1/2 that of Model 2: ( 
can stay high enough for the MRI to operate until very 
near the midplane, where Ohmic diffusion begins to dom- 
inate. An open question is whether or not a dead slice 
would be present in a fully 3-D simulation with identical 
vertical profiles of p, ( and /3 to our disk: the thickness 
of the dead slice is smaller than the MRI wavelength by 
a factor of 2-10, so turbulence would very likely erase it. 

We have seen that MRI-active disks may have com- 
plex vertical structure, with multiple layers of turbulent 
and non- turbulent zones. In the next section we explore 
the overall mass fiow through the disk and discuss its 
evolution on million-year timescales. 



4.2. Mass Flow 

As might be expected for a layered disk with a dead 
zone, radial mass transport is not in steady state for 
either Model 1 or Model 2. The left-hand panel of Figure 
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Fig. 3. — Plots of viscosity u (top), turbulent stress coefficient a (middle), and ionization rate C, (bottom) for the inner 4 AU of Model 
2 at 1 Myr. In the inner 3 AU, a stable corona where ambipolar diffusion shuts down MRI sits on top of the active layer. Note the thinness 
of the active layer over the pileup at 1 AU. The dead slice seen in Figure^ is also visible between 1.2 AU and 1.7 AU. 
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m shows 



M{R) = 2ttRY.vr 



(54) 



for both models after 1 Myr of evolution, where vr is a 
density-weighted, vertically averaged gas radial velocity. 
The convention in Equation[5^is that vr is positive when 
gas flows toward the star and negative when gas flows 
away from the star. In the inner part of the disk, where 
gas flows inward, the highest accretion rates occur (a) at 
the outer edge of the dead zone (at 21 AU for Model 1 and 
16 AU for Model 2), and (b) at the inner disk boundary. 
There is a clear drop in M{R) associated with the dead 
zone. (The noisy M profile where the disk has two or 
more layers is a reflection of the root-finding tolerance 
in our Newton-Raphson algorithm. Fluctuations in the 
location of the boundary between dead and active zones 
are random and average out over time.) 

The potential for the dead zone to become gravitation- 
ally unstabl e, given enough time to accumulate mass, is 



clear ( jGammie 1996 Armitage et al. 2001 Zhu et al. 
|2010[ )— though we see only a slow, steady density growth 
at ^ 1 AU-4 AU over the course of 3 Myr in Model 
2 (Figure |5]). We will examine the potential for gravi- 
tational instability in high-mass, MRI-active disks in a 
forthcoming paper. One important caveat is that the 
accumulation of mass in the dead zone is unstable to 



the Rossby wave i nstability ( Hawley|1987[|Li et al.|2001[ 
Meheut et al.|[2Qr2J , which triggers spiral density waves. 
Without aziniuthal information in our model, we can- 
not model the effect of the Rossby wave instability on 
dead-zone overdensity. Ultimately, the overdensity may 
not survive and m ay break up into la rge-scale vortices 
(Papaloizou & Pr ingle 1985 Lovelace et al. 1999). 

At t = 1 Myr, there is turn-around in the accretion flow 
near 40 AU in each model. The outward mass flow out- 
side 40 AU is necessary for overall conservation of angular 
momentum as material in the inner disk moves toward 
the star. The fact that there must be a change in the ac- 
cretion flow direction also follows directly from the diffu- 
sive n ature of T-Tauri disks (e.g. Lynden-Bell & Pringle 



1974), which lack a circumstellar envelope to feed steady 



state inward accretion. In both models, the disk has ex- 
panded beyond its initial 70 AU size after 1 Myr. Note 
that the 40 AU region of the disk does not drain, as the 
"turn-around radius" where the accretion flow changes 
direction moves inward with time (right-hand panel of 
Figure [4| . We heartily discourage the use of a steady- 
state accretion rate for any reasonable description of an- 
gular momentum transport in T-Tauri disks. However, 
for predicting observables to order-of-magnitude using a 
static, non-time-varying disk model, we suggest approxi- 
mating the accretional heating in MRI-active disks with 
a constant M ~ 10~^Mq yr~^. Figure [i] suggests that at 
a given time, M has a modest dependence on disk mass. 

How consistent are our modeled accretion rates with 
observations? An oft-quoted value of M^, the T-Tauri 
star accretion rate, is 10~^Mq yT~^ ( Sicilia- Aguilar et al 
[2004 ). Accretion rates in transitional disks are roughly 
3xlO"^M0 yr-^ ( [Espaillat et al.|2012[ ). Our simulations 
achieve an accretion rate of roughly ~ 10~^Mq yr~^, 
more consistent with the median transitional disk accre- 
tion rate. Reproducing high accretion rates in the inner 
regions of MRI-active disks is a delicate balance between 



magnetic field strength. X-ray ionization, grain size, and 
gas/small grain mass ratio. For example, we performed 
simulations using the same disk masses and X-ray lumi- 
nosity (O.O15M0 and 0.03 M©; 2 x 10^° erg s-\ respec- 
tively) but a grain size of 0.1/im instead of 1/im and a 
standard gas/small grain mass ratio of 100, and quenched 
all M RI activity in the di sk entirely. However, as mod- 
els by Zsom et al.| (2011) show grain growth and set- 
tling within 1000 years, it is reasonable to assume the 
gas/small grain mass ratio has evolved from the inter- 
stellar value by the T-Tauri phase. In our simulations, a 
stronger magnetic field (lower value of P at the midplane) 
also suppressed MRI-driven accretion, which occurs in re- 
gions of the disk dominated by ambipolar diffusion. Most 
simulations of non-ideal MRI-driven accretion ha ve diffi- 
culty reaching the 10~^ Mq yr~^ benchmark (e.g. |Zhu et| 
al.|2010[[Bai|20 11), though there IS a substantial amount 
of sca tter in both T-Tauri and transitional disk accretion 
rates ( Romero et al.|[2QT2 ). 



4.3. Simple prescription for accretional heating in an 
MRI-active disk 

Previous studies have often relied on static, non- 
evolving disk models to connect measured line fluxes, ve- 
lociti es or spectral indices to physical propertie s of disks 



(e.g. D^Alessio et al. 2006 Pinte et al. 2010). Unfor- 



tunately the constant-a prescription lor turbulent angu- 
lar momentum transport leads to unphysical thermody- 
namic descriptions of disk midplanes. Predicted mid- 
plane temperatures in a constant-a model, where most 
turbulence is concentrated at the midplane, are too high. 
The poor match of the constant-a model to realistic pro- 
tostellar disk accretion is definitely a problem for obser- 
vations that trace disk midplanes, such as sub-millimeter 
measurements of continuum emission from large grains 
or surveys of rare molecules like HCO+ or HD. How- 
ever, the problem also affects observations that trace 
surface layers, such as Spitzer emission lines, infrared 
SEDs and sub-millimeter maps of abundant gases like 
CO. If a significant subset of T-Tauri disks rely on MRI 
to drive accretion, overestimating their midplane tem- 
peratures leads to overpredicting the photosphere and 
pressure scale heights, underpredicting the optical depths 
of low-excitation lines such as CO (J = 2 ^ 1), and 
possibly underpredicting the rates of grain settling and 
growth. 

Here we use our simulation results to present a simple 
prescription for a in an MRI-active disk. First we deter- 
mine the depth of the MRI-active layer. Figure [6] shows 
the surface density of the active layer as a function of ra- 
dius for four time snapshots of Model 1 (left) and Model 
2 (right). Clearly, the functional form of i]active(^) is 
the same for both models and does not vary significantly 
with time. Near the star, where X-ray ionization is im- 
portant, ^active IS high but falls off quickly as X-ray irra- 
diation declines (see Equation 29). Over the dead zone. 



in the region where only cosmic-ray ionization is impor- 
tant, i;active(^) iucreascs as the magnetic field weakens, 
shrinking the inactive corona. I]active(^) reaches its max- 
imum at the outer edge of the dead zone. Here, where 
the entire vertical column is active, the decreasing depth 
of the active layer simply reflects the fact that S(i?) is a 
decreasing function. 
Although 

^active docs Vary with R in the parts of the 
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1 Myr 




R (AU) Time (yr) 

Fig. 4. — Left: A plot of M{R) for both disk models after 1 Myr of disk evolution reveals that low-mass MRI-driven accretion disks have 
not reached steady state after 1 Myr. As required by conservation of angular momentum, there is a turn-around in the mass flow between 
inward accretion and outward decretion. The highest accretion rates in the inner disk, where mass moves toward the star, occur on either 
side of the dead zone. Right: The location of the turn-around radius as a function of time for both models. The turn-around radius moves 
inward over time. 
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Fig. 5. — Volume density as a function of {R, z) at two timepoints during disk evolution. In the more massive Model 2, the density in 
the pileup at 1 AU noticeably grows with time. Top left: Model 1 (disk mass O.OISM©), 10"^ years of simulation time (star age 0.11 Myr). 
Top right: Model 1 after 1 Myr of simulation time. Bottom left: Model 2 after 10^ years of simulation time. Bottom right: Model 2 
after 1 Myr of simulation time. 



disk with a dead zone, the variation is less than a factor of 
two for both models presented here. Similarly, Sactivel^) 
is slightly higher for Model 2 than for Model 1, suggest- 
ing that the active column may increase modestly with 
disk mass, ^active ^ 10 g cm~^ is a good approximation 
for low-mass disks with surface density less than about 
four times the MMSN. The semi-analytical viscosity pre- 
scription is then simple: 

1. Set the depth of the active column to Sactive ^ 
10 g cm~^. Set a ~ 0.01 in the active column and 
a ~ 10~^ in the dead zone. While there may be a 
corona on the inner disk surface, it contains very 



little mass and may safely be ignored as long as the 
disk model includes stellar heating. 

2. Smooth the transition between the active layer and 
the dead zone if desired. 

3. The entire vertical column will be active where S ~ 
20 g cm~^, such that the two active layers meet at 
the midplane. Outside the outer radius of the dead 
zone, where S ^ 20 g cm~^, use a ~ 0.01. 

4. Calculate viscosity u at each (i?, z) in the disk using 
Equation [l] 
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Fig. 6. — Eactive(^) for four time snapshots of Model 1 (left) and Model 2 (right). The active column depth as a function of radius 
changes very little as the disk evolves, making it possible to approximate Eactive(^) for use in a semi-analytical viscosity prescription. 



5. THERMAL EVOLUTION OF MRLACTIVE DISKS 

Now that we understand how mass moves through a 
disk where accretion is driven by MRI, we turn our atten- 
tion to the thermal structure and evolution of the disk. 
Here we see some important differences from disk models 
tha t use a constant-a viscosity prescription. In Section 
|5.1| we discuss how different parts of MRI- active disks 
heat up and coo l off over time (Question 4 in Introduc- 
tion). In Section j5.2| we investigate the location of the ice 
line and its motion through an MRI- active disk (Ques- 
tion 5 in Introduction). 

5.1. Disk Heating and Cooling 

The top panels of Figure [7| show the surface and mid- 
plane temperatures of Models 1 and 2 after 10 thousand 
years of evolution (left) and 1 Myr of evolution (right). 
First, note the obvious feature that the surface is far 
hotter than the midplane throughout most of the disk. 
Disk models with constant-a viscosity prescriptions usu- 
ally have warm midplanes, T > 100 K, in the inner 5 AU. 
In the inner 1 AU of constant-a models, the midplane 
temperature can even exceed the surface temperature, 
approaching 1000 K (e.g. j lersant et al. 2001; D'Alessio] 
et al.|2006 Dodson- Robinson et al. 2009). In our models, 
there is so little turbulent energy generated in the dead 
zone that the disk midplane falls to 20 K — the ambient 
temperature of the remnant molecular cloud surround- 
ing the disk. Here we assume the disk is optically thin 
to long- wavelength radiation from the ambient cloud and 
cannot cool below the ambient temperature. Only in the 
inner ~ 2 AU does residual mass transfer by large-scale 
magnetic fields lift the midplane temperature above the 
minimum value. 

Moving outward through the disk, the midplane tem- 
perature rises modestly until it equalizes with the surface 
temperature where the disk becomes optically thin to 
stellar radiation. Both of our models, despite their mod- 
est masses (O.OISMq and O.OSM©), are optically thick 
to stellar radiation inside ^ 70 AU, decoupling the sur- 
face and mi dplane temperatures . Passively heated disk 
models (e.g. Woitke et al."2009) are therefore poor ap- 
proximations to the midplane temperatures of our model 
disks, as are const ant-a models in which T decreases with 
R at the midplane. Instea d, co upling the simple viscos- 
ity prescription in Section 4.3 with a radiative transfer 
scheme is preferable for modeling observables that trace 
the inner ~ 70 AU. 



Figure [T] shows that the disk surface cools off over time, 
as expected for any disk being irradiated by a T-Tauri 
star evolving along the Hayashi track. The surface cool- 
ing affects the disk structure as follows: 

1. The overall photosphere height of the disk de- 
creases with time (Equation 39 Figures l2| [s] and 

III); 

2. The viscosity in the surface layers decreases with 
time (Equation 31 Figure |2|. 



The decrease in surface viscosity with time is due to the 
increasing density in the surface layers as the disk cools. 

Despite the fact that the disk surface cools with time, 
the tendency of the MRI to pile up mass unevenly leads 
the temperature to increase with time in certain parts 
of the disk. The bottom panels of Figure [7| show the 
temperature increase in the pileup at 1 AU in Model 
2. Though modest, the temperature increase does affect 
the location of the ice line, which we discuss in the next 
section. 

5.2. The Ice Line in MRI- Active Disks 

Ice forms in protoplanetary disks when the tempera- 
ture falls below 145-17 K, depending on the water va- 
por's par tial pressure (Podolak fc Zucker||2004[ Lecar et 
al. 2006). Observations of the outer asteroid belt place 
the ice Ime in today's Solar System at 2.7 AU. Theoret- 
ical estimates of the ice line location in the solar nebula 



place it a minimum of 0.6 AU from the sun (|Davis|[2005 
Garaud fc Lin||20Q7 ) and a maximum of 6 AU at the be- 
ginning of the 1 -'I'auri phase, moving in ward as the disk 
evolves ( Dodson- Robinson et al. 2009). In constant-a 
disks, the ice line moves inward with time as the opti- 
cally thick disk radiates away its accretion energy. At 
late times, however, the inner disk loses enough mass 
to become optically thin to stellar irradiat ion, causing 
the ice line to move o utward with time (Garaud fc Lin 



Oka et al. 2011). The more massive a constant-ce 



h e high e r"^ its midplane temperature w ill b e. Tak 



2007 

disk^ . 

ing the^ ^ ( |2005[ ), [Garaud fc Liii| (|2007l ) and jOka et 
al. (2011) disk models up to a reasonable planet-iormmg 
mass of O.O4M0 or higher (Thommes et al. 2008) would 
push their ice lines outside the terrestrial plane t^orming 
region, more in agree ment with the results of Dodson- 
Robinson et"aL] ( |2009| ). 

In an MKI- active disk, the overall thermal evolution 
is determined by the dimming of the parent star, which 
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Fig. 7. — Although the surface of the disk cools over time, as the T-Tauri star moves down the Hayashi trach, parts of the midplane 
of an MRI-active disk may heat up over time as mass piles up in the dead zone. Top left: Temperature at the surface (top curves) 
and the midplane (bottom curves) of Model 1 (red dotted) and Model 2 (black solid) after 10,000 yr. Top right: Surface and midplane 
temperatures of Models 1 and 2 after 1 Myr. Note the increase in midplane temperature at the edge of the dead zone in each model. 
Bottom left: Contour plot of T(R, z) in the inner 4 AU of Model 2 after 10,000 yr of evolution. Temperature units are Kelvins. Bottom 
right: Contour plot of T(R, z) in the inner 4 AU of Model 2 after 1 Myr of evolution. Note how the pileup of material centered at 1 AU 
(see Figure |5| has heated up over time. 



causes viscosity in the activ e surface layers to decrease 
with time (see Section 5.1). Figure [s] shows the 2-D 
structure of the ice hue in Model 2 at timepoints 100 yr, 
1000 yr, 10,000 yr, 100,000 yr, 1 Myr and 3 Myr. At early 
times, we re cover the ^^tw o-branch" structure of the ice 
line seen by Davis (2005), in which a nearly horizontal 
ice line divides the hot surface from the cool interior and 
a midplane ice line separates the warm midplane near 
the star from the cool midplane far away from the star. 
After 10,000 yr of evolution, the inner edge of the disk 
at 0.5 AU cools enough for ice to freeze at the midplane, 
causing a "pinch-off" in the midplane branch of the ice 
line that leaves an H2O gas bubble at 0.6 AU. This pinch- 
off is probably a boundary effect: since mass flows freely 
from our inner boundary at 0.5 AU onto the star, the disk 
near our inner boundary loses mass and cools quickly. 

By 1 Myr, mass loss from the inner edge of our grid 
combined with less activity in the surface layers cools the 
disk enough to push the midplane ice line inside 0.5 AU, 
the inner boundary of our computation. Here our results 
are st il l consi stent wit h the low- mass models of iDavisI 
( [20051 ) ^ [Gar^d fc Linj ( [2QQ7| ) and |Oka et aI| ( [2QTT] ), who 
all show the ice line movmg inside 1 AU for accretion 
rates M ^ 10~^^Mq yr~^. An impor t ant d ifference be- 
twee n our mode ls and thos e of Davis (|2005|), [Garaud & 
Lin (2007) and Oka et al. ( 20lip is that the accretion 
rate does not have to drop extremely low for the terres- 
trial planet-forming region to cool enough to freeze ice: 



both Model 1 and Model 2 have active layers that drive 
M ~ 1O~^M0 yr~^ through the disk surface. Our model 
disks pr edict colder midp lanes in the inner 10 AU than 



those of Terquem 
assumed here (a 

or ajj = 10~^ in Figure 3 of Terquem 



(|2008|) due to the lower minimum a 

3 



10 for our work vs. ap = 10 



(|2Q08|), wh ere a 



the value of a in the dead zone). Unlike|TSquem (2008[), 
our disk has a colder midplane than surface because we 
include the effects of stellar heating. 

The MRI can also create local regions that heat up 
with time not because the disk becomes optically thin, 
but because mass piles up in the dead zone. Between 
1 Myr and 3 Myr of evolution, the pileup at 1 AU of 
Model 2 heats up enough for the midplane ice line to 
reappear — a real physical effect. Model 1, which does 
not show any such pileup, follows a similar evolutionary 
path as Model 2 up to 1 Myr. In an MRI-active disk, 
water ice in the terrestrial planet-forming region may be 
transient. Understanding which parts of the inner disk 
have water ice available for planet formation requires a 
careful comparison of the planet esimal growth timescale, 
the star cooling timescale and the growth timescale of 
any pileups deposited in the dead zone. In Model 2, the 
pileup grows and the star dims on a timescale similar to 
the disk lifetime. We will examine high-mass MRI-active 
disks, in which dead-zone pileups can grow much more 
quickly, in a forthcoming paper. 
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Fig. 8. — In the inner 2 AU of Model 2, we see that the disk surface coohng and midplane heating both affect the structure of the ice hue. 
Each plot shows T(R,z) with the ice line at 160 K indicated by the blue contour (temperature units are Kelvins). Simulation times are 
100 yr (top left), 1000 yr (middle left), 10,000 yr (bottom left), 100,000 yr (top right), 1 Myr (middle right) and 3 Myr (bottom 
right). Recall that the star age is already 0.1 Myr at the beginning of the simulation. 



6. MODEL LIMITATIONS 

Although our model incorporates much of the physics 
of disk evolution during the T-Tauri phases of stars, our 
long-timescale computation requires a number of sim- 
plifications. Our disk accretion model suffers from the 
following limitations, which may affect our conclusions: 



1. 



We model our disks as 1+ld instead of 3d. We as- 
sume radial symmetry, and that the vertical struc- 
ture is not coupled to the radial mass transport. 
This ignores the possibility of accretion driven by 
non-axisymmetric instabilities not represented by 
our parameterization of stress from MH D turbu- 
lence , such as the Rossby wave instabilit y (Hawley 
1987[|Li et a I]|2QQl| peheut et al.|[2QT2l ). We may 



then underestimate accretion rates m some regions 
of the disk. 

2. We do not model the effects of gravitational insta- 
bility on momen tum transport, as is done by |Mar-| 



tin et al. (12012b), but instead limit our simulations 
to those disks which are gravitationally stable for 
their entire lifetime. To test whether a disk is gravi- 
tationally stable to axisymmetric perturbations, we 



calculate the Toomre Q parameter at every radial 
grid point. Stability requires that 



3. 



> Qc 



1. 



(55) 



Note that for disks more massive than those we 
have simulated, the disk does eventually become 
gravitationally unstable to axisymmetri c pe rturba- 
tions within the dead zone (see Section 4.2). 



We do not model the effects of Hall diffusivity, 
which is the least well understood magnetic diffu- 
sion regime. However, it is unlikely to alter either 
the conditions required for MRI or the strength 
of the t urbulence where O hmic dissipation is also 
present ( |Sano fc Stone|[2QQ2p . 



4. We make several assumptions about the magnetic 
field strength and presence of grains. The activity 
of MRI driven turbulence is very sensitive to these 
parameters. First, we deviate from the standard 
interstellar gas/small grain mass ratio of 100 and 
instead assume a ratio of 1000. This is equivalent 
to assuming that 90% of the grain mass is either 
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in grains larger t han Ijam (Oliveira et al. 2010 
Scaife et al. 2012), which do not signiticantly at- 
tect e lectron density, or has settled below the dead 



zone. [Mohanty et al.| ( |2013[ ) find that grain deple 
tion through growth or settling is required to ac- 
count for the observed accretion rates of low mass 
protostars. Without the high gas/small grain mass 
ratio, the inner Ohmic dead zone and the upper 
Ambipolar dead zone overlap, producing a passive 
thermal structure for much of the radial extent of 



the disk (see Section 4.1). 



We also assume a vertically constant magnetic 
pressure with a midplane plasma /3 of 1000. This 
is at the upper end of the range for m idplane 
plasma (3 found by [Fromang fc Nelsonj ( |2006[ ) 
in global MHD simulations with saturated turbu- 
lence. Without a sufficiently weak magnetic field, 
the ambipolar diffusivity can be large enough that 
there is once again an overlap in the dead zones 
and an at least partially passive disk. 

5. We assume a zero stress boundary c ondi tion at the 
inner radius of our disks (see Section 3.3 ). This cre- 
ates a non-realistic boundary effect m which mass 
is rapidly depleted from the inner annuli. One re- 
sult of this can be seen in Figure [Sj as the midplane 
within 0.6 AU becomes cold due to the loss of sur- 
face density. 

Due to computational convergence difficulties in our 
model, at any particular timestep there are a few annuli 
with vertical structures which do not meet one or more 



of our midplane boundary conditions (Equations 50 , 51 



and 



521), due to the discontinuous nature of the MKl anc 
We do not consider this to be a limitation on 



opacity 

our simulations, as the unsolved annuli's contribution to 
the viscosity profile is smoothed with a Gaussian filter 
before being used to update the surface density profile 
(Equation 46). Although there are always a few "bad" 
annuli present, at any particular radius the lack of con- 
vergence for the vertical structure persists for only a few 
timesteps. The unconverged annuli can be seen as oc- 
casional mcongruous vertical bars in our contour plots 
(Figures 



mcongruous verti 
|2l[3l[5|[7land[8|. 



7. CONCLUSIONS 

In the Introduction, we asked five questions about the 
structure and evolution of MRI-active disks. Here we 
summarize our findings and answer each question: 

1. How do the relative sizes of the dead zone and ac- 
tive layers change over time? 

The radial size of the dead zone is almost con- 
stant in time, while the vertical height of the dead 
zone shrinks over time. What was surprising about 
our results was not the evolution of the dead zone, 
but the complexity of the disk structure. Between 
1.2 AU and 1.7 AU, Model 2 has a five-layer struc- 
ture throughout most of its evolution: inactive 
corona, active layer, d ead slice, active layer, dead 
midplane (see Section 4.1 and Figures fl] and p|. 
Note, however, that detailed 3-D simulations would 
likely show no dead slice since the MRI wavelength 



in the active layers bracketing the dead slice is of or- 
der the dead slice thickness. The lower-mass Model 
1 has at most three layers: inactive corona, active 
layer, dead midplane. Throughout this work, we 
have seen that increasing disk mass leads to in- 
creasing complexity in the disk structure and ac- 
cretion fiow. 

Finally, the radial size of the dead zone was 
somewhat higher than predicted in previous work: 
16 AU for Model 1 and 21 AU for Model 2. Previ- 
ous papers reporting a ^ 5 AU dead zone used the 
MMSN (e g. | M atsumu ra fc Pudritz|2003[ [Salmeron 
fc Wardle|2008( rFlaig et al .|2012[) , but the modestly 
higher surface dehsities of our disks expanded the 
dead zone. There is some evidence that the solar 
nebula had a large dead zone consistent with our 
findings — the giant planets have an atmospheric 
composition gradient that, if primordial, would 
have diffused on million- ye ar timescales if not pro - 
tected by a dead zone (Nelson fc Gressel 2010). 
One caveat, though, is that we have assumed that 
the stellar X-ray ffux is constant in time, as is the 
cosmic ray ffux. A decreasing stellar X-ray ffux, 
which ionizes mainly the inner ~ 3 AU of the disk 
surface, might erase the dead slice over time, while 
changing the cosmic ray ffux as the ambient molec- 
ular cloud disperses would certainly change the ra- 
dial extent of the dead zone. 

2. How does M vary with radius and time? 

Throughout the disk evolution, \M\ is highest at 
the inner and outer boundaries of the dead zone. 
The lowest \M\ in the inner disk, where gas ffows 
toward the star, is in the middle of the dead zone. 
|M| is about 50% higher for Model 2 than Model 1, 
suggesting that higher-mass disks support higher 
MRI-driven accretion rates — though the increase 
in M with Mdisk is modest. \M{t)\ decreases ex- 
tremely slowly: though the depth of the active layer 
does not change with time (Figure [6|, the viscosity 
in the active layer drops modestly as the star cools. 

As required by conservation of angular momentum, 
both model disks expand as they evolve, creating a 
turn-around in the accretion ffow. The turn-around 
radius is almost entirely determined by Rout at 
t = and moves inward as the disk evolves. Here, 
with Rout = 70 AU at t = 0, the turn- around radius 
eventually reaches 40 AU after 1 Myr of evolution. 
Note that the turn-around radius moves steadily 
inward in both models and does not converge to- 
ward a particular location (Figure [4| . One expects 
the turn-around radius to move steadily inward be- 
cause mass must continually join the outward "de- 
cretion" ffow in order to keep transporting angular 
momentum outward. 

Our models predict M ^ IO'^Mq yr ^ in the 
planet-forming region of the disk, but the exact 
value of M depends on many free parameters such 
as grain size, gas/small grain mass ratio and mag- 
netic ffeld strength. We have not attempted an 
exhaustive parameter study of M as a function of 
all variables. We merely note that for a gas/small 
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grain mass ratio of 100 and a grain size of 0.1/im, 
all MRI activity in the disk was suppressed. Like- 
wise, decreasing plasma /3 at the midplane from 
1000 to 100 suppressed MRI turbulence, though 
not as severely as small grains. 

3. How can disk modelers parameterize heating in the 
active layers and dead zone without resorting to a 
3-D MHD simulation? 

The value of non-evolving, "snapshot" disk mod- 
els is inargua ble, p articularly for modeling observ- 
ables. Section |4]3] presents a simple modification of 
the standard, constant-a irradiated disk that ap- 
proximates the thermal structure of the disk where 
the dead zone is present. Simply set the surface 
density of the active layer to 10_g cm~^ and use 
a ^ 0.01 (see Figures [l] |3] and [6|. For the dead 
zone, use a ~ 10~^. For numerical models, we 
recommend smoothing the transition between the 
dead zone and the active layer. While the depth of 
the active layer does vary across the dead zone, its 
variation is at most a factor of two in a given disk. 
The fact that the active layer depth is almost con- 
stant in time for both Model 1 and Model 2 makes 
our simple prescription applicable to any stage of 
T-Tauri disk evolution, provided that stellar irra- 
diation is included. 

4. Does the disk midplane heat up or cool off with 
time? 

In both models, the midplane temperature varies 
little with time, but the disk surface cools as the 
star evolves down the Hayashi track. Despite the 
modest masses of our model disks, the optical 
depth of both to stellar irradiation is enough to 
thermally decouple the surface from the midplane. 
Most of the dead zone is so lacking in energy gener- 
ation that it falls to the assumed ambient tempera- 
ture of surrounding molecular cloud material, 20 K 
in these models. Since the radial extent of the dead 
zone changes little with time, the disk midplane 
temperature remains static except at > 60 AU, 
where the disk thins enough over time to become 
optically thin to stellar irradiation (Figure [t]). 

In the midplane, there are two possible locations 
where the temperature is not static but increases 
with time (Figure [t]). The first is at the outer 
edge of the dead zone. A slight decrease in surface 
density with time pushes the dead zone boundary 
modestly inward, allowing material at the edge of 
the dead zone to become turbulent and heat up. 
The other location of increasing temperature with 
time is the pileup at 1 AU of Model 2. Lower- 
mass Model 1 does not develop any such pileups on 
million- year timescales. To the extent that MRI- 



driven accretion deposits piles of material in the 
dead zone, the disk midplane may heat up. 

Note, however, that the thermal properties of the 
disk depend on the ionizing radiation it receives. 
Near the star, the disk structure would evolve if 
the X-ray fiux were to change with time. In future 
work, it would be interesting to let Lx scale with 
bolometric luminosity on the Hayashi track. A 
spike in cosmic ray ionization from a near by super- 
nova woul d affect the global disk structure ( Fatuzzo 
et al.|[2QQ6 ). 



5. Where is the ice line in an MRI- active disk and how 
does its location change over time? 

Due to the paucity of energy generation in the dead 
zone, the midplane ice line falls somewhere inside 
the terrestrial planet-forming region. In our mod- 
els, the midplane ice line actually moves off the in- 
ner edge of our grid at 0.5 AU after 10^ years of disk 
evolution, reappearing in Model 2 after the pileup 
at 1 AU reheats the midplane. Despite the dif- 
ferent physics used in computing the ice lin e loca- 
tion, our results roughly agr ee wi th those of Davis 
( [20051 ) ^ [Garaud fc Lin| ( [2QQ7| ) and jOkaet al.||2QTT[ 
in that the midplane ice line is inside 1 AU tor most 
of the disk's evolution. Determining whether and 
when ice is available for terrestrial planet forma- 
tion requires carefully comparing the disk cooling 
timescale, planetesimal growth timescale and the 
timescale on which MRI-deposited pileups grow. 
An important difference between our model and 
standard, constant-a disk models is that the ac- 
cretion rate does not have to drop extremely low 
to move the ice line inside 1 AU: throughout their 
evolution, our disks have M > lO^^M© yr~^ mov- 
ing through the active layers that sandwich the icy 
inner disk. 

Here we have presented the first analysis of the struc- 
ture and evolution of an entire MRI-active disk on 
million-year timescales. While we have chosen to focus 
on low-mass disks in this work in order to compare with 
previous studies, we will expand our analysis to include 
high-mass, planet-forming disks in a forthcoming paper. 
Already, at Mdisk = O.O3M0 — ^just short of the minimum 
0.04M(7) req uired for giant planet formation (Thommes 
et al. 2008") — Model 2 exhibits some new features such 
as the split active layer and the re-heating midplane. 
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TABLE 1 

Parameters in protostellar disk model 



Name 



Value 



Description 



Reference 



Ionization parameters 



/■surf 
^CR 


10-17 s-1 


Cosmic ray ionization rate at disk surface 


^CR 


96 g cm-^ 


Cosmic ray penetration depth 


Lx 


2 X 10^0 erg s-^ 


Stellar X-ray luminosity 


Ci 


6.0 X 10-12 s-i 


Ionization rate coefficient for absorbed X-rays 


C2 


1.0 X 10-1^ s-i 


Ionization rate coefficient for scattered X-rays 


Ai 


2.5 X 10-3 g cm-2 


Penetration depth of absorbed X-rays 


A2 


1.2 g cm-2 


Penetration depth of scattered X-rays 


PI 


0.4 


Exponent of absorbed X-ray attenuation 


P2 


0.65 


Exponent of scattered X-ray attenuation 


Ambient medium 






20 K 


Ambient temperature set by remnant molecular cloud 


Dust ^ 


grains 




pgr 


3.0 g cm-3 


Internal grain density 


G/S 


1000 


Gas/small grain mass ratio 


a 


1 /im 


Grain size 


Gas composition 




Nmg 




Magnesium abundance in disk gas 


LL 


2.33 g mol-i 


Mean molar weight 


Maxwell stresses 




10-^ 


Minimum stress in dead zone from large-scale fields 



Umebayashi &: Nakano ("iQSlt 
Umebayashi &; JNakano (lllHTl 
^armire et al. (2000) 
Kai &; Goodman (20091 
Hai &; Goodman 
Bai &; Goodman 
Hai &; Goodman 
Bai <fe Goodman 



(2001) 
(2001) 
(200^) 
(2UU^) 



Hai Goodman] ( pOUU ) 



[Peretto etal] ( |2010| ) 

standard 

augmented from standard 100 to approximate grain growth 
|01iveira et al.| ( |2010| ) 



Turner & Sano 
standard 



( 2008| ) 



Turner Sano| ( |2008| ) 



